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Nonlinear transport in the one dimensional Hubbard model at half-filling under a finite bias 
voltage is investigated by the adaptive time-dependent density matrix renormalization group 
method. For repulsive on-site interaction, dielectric breakdown of the Mott insulating ground 
state to a current-carrying nonequilibrium steady state is clearly observed when the voltage 
exceeds the charge gap. It is found that by increasing the voltage further the current-voltage 
characteristics are scaled only by the charge gap and the scaling curve exhibits almost linear 
dependence on the voltage whose slope is suppressed by the electron correlation. In the case 
of attractive interaction the linear conductance is the perfect one 2e^ /h which agrees with the 
prediction by the Luttinger liquid theory. 

KEYWORDS: Hubbard model, Mott insulator, dielectric breakdown, nonequilibrium steady state, adaptive 
time-dependent DMRG 



Emergence of rich variety of different states of mat- 
ters is a consequence of electron-electron interaction. 
The Hubbard model is a prototypical interacting electron 
system and shows different phases depending on lattice 
structure, filling and interaction. When the band is half- 
filled and the Coulomb repulsion is sufficiently strong, 
the charge excitations involve a finite energy gap Ac, 
and this fact is a manifestation of the Mott insulating 
ground state. In one dimension (1-D) intriguing proper- 
ties of the model including the Mott transition have been 
clarified by various analytic approaches: the Tomonaga- 
Luttinger liquid theory, the Bethe Ansatz and the con- 
formal field theory.^' Therefore concerning equilibrium 
properties one can say that the 1-D Hubbard model is 
the best studied model in depth. 

Instead of chemical doping which is commonly used 
to realize metal-insulator transitions, one can also ap- 
ply a bias voltage to break an insulating phase. ^^ This 
process, dielectric breakdown of a Mott insulator, may 
be called as a nonequilibrium metal-insulator transition. 
However, no systematic theoretical study on the break- 
down of Mott insulators have been performed due to the 
difficulty to treat the nonequilibrium states of strongly 
correlated systems. 

Recently the adaptive time-dependent density matrix 
renormalization group (TdDMRG) algorithm was devel- 
oped,^' which is an extension of the DMRC^' method 
to time-dependent problems. This technique has been 
used as a powerful numerical approach to nonequilibrium 
problems in one spatial dimension with strong correla- 
tion, such as single quantum dot system under finite bias 
voltages^' and the interacting resonant level model. ^^ 

Oka and Aoki utilized the TdDMRG method to study 
the breakdown of the Mott insulating phase of the 1- 
D Hubbard model driven by an external electric field. ^' 
They demonstrated that the phenomenological expres- 
sion for the transition probability which is obtained by 
replacing the band gap in the Landau-Zener formula by 

* E-mail address: kirino@issp.u-tokyo.ac.jp 



the many-body charge gap Ac is consistent with the U 
dependence of the threshold electric field. They could dis- 
cuss, however, only the threshold and it is necessary to in- 
vestigate current- voltage (I-V) characteristics beyond the 
breakdown to elucidate nature of nonequilibrium steady 
states of the 1-D Hubbard model. 

In this Letter the 1-D Hubbard model with a finite 
bias voltage is studied by the TdDMRG method which 
enables us to obtain for the first time reliable numerical 
results on currents. These results are clear manifesta- 
tion that various nonequilibrium phenomena in strongly 
correlated 1-D systems have become within the reach of 
theoretical investigations. 

We determine the I-V characteristics for the repulsive 
1-D Hubbard model at half-filling, and show that nonzero 
steady current appears when the bias voltage exceeds Ac. 
By increasing the voltage beyond Ac current is scaled 
only by the Ac if the voltage is small compared to the 
band width and the scaling curve has almost linear re- 
gion with its slope suppressed by the correlation effect 
compared to a band insulator. Concerning the attractive 
case, low energy properties of the 1-D Hubbard model 
at half-filling are classified into the Luther-Emery liq- 
uids which are characterized by gapless charge excita- 
tions and gapful spin excitations, and the ground state 
of the model has two degenerated quasi long range or- 
ders, superconducting and CDW ones.*-' We show that 
the linear conductance of the attractive Hubbard model 
is precisely given by the perfect conductance 2e^//i. 

We consider the 1-D Hubbard chain at half-filling with 
an applied DC voltage. Our main target in this Letter is 
the nonequilibrium steady states of the system at T = 0. 
Although the electric potential inside the system should 
be determined self-consistently, we neglect the change of 
the electric potential due to the charge redistribution. 
Since the TdDMRG method is applicable only to finite 
systems, it is important to reduce system size dependence 
of the results. Thus we take for the voltage term the sim- 
plest model in which the potential difference is confined 
solely to the central bond. In doing so the current- voltage 
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characteristics can be addressed by the finite system cal- 
culation as we will show below. 

In order to realize nonequilibrium steady states in the 
numerical calculation, we first obtain the ground state 
wave function of the system without the voltage by the 
standard DMRG method. Then we calculate the time 
evolution of the wave function after the switching-on of 
the bias voltage using the TdDMRG algorithm.^' The 
nonequilibrium steady state is described by the wave 
function after some transient period. Putting the above 
information, the Hamiltonian is written as 



H{t) = Hl + Hr - t' ^(4c,, + h.c.) 
Ha = -t ^ (cLcj+iCT + h.c.) 



(1) 



+ ^Y. 4t4lC^iC,^ (« ^L,R), (2) 

where L{R) represents the left (right) half of the system, 
Cio- annihilates an electron with spin a at ith site, r is the 
time variable, t the hopping amplitude, U the Coulomb 
energy, V the applied voltage and Na = J2tea,a4a^i'^- 
/(r) is the index of rightmost (leftmost) site in L{R) 
and t' is the hopping between the Ith and rth sites. In 
this Letter we concentrate on i' = i case for simplicity. 
The bias voltage is turned on according to the smoothed 
step function 9{t) = (1 + cxp [(tq — t)/ti])~^ in order 
to mimic adiabatic switching-on and to soften transient 
behaviors. We fix tq = 4fi,/t and ri = h/t throughout 
this Letter. 

At half-filling an electron-hole transformation for one 
species of spin, known as the Shiba transformation,^' 



Cjt -^ cjt, cji -^{-lyc 



w 



(3) 



maps the U > {U < 0) Hubbard model to the U < 
{U > 0) one. The charge and spin degrees of freedom are 
interchanged by this transformation. For example, the 
bias voltage term is transformed to the Zeeman term, 
and the current operator is replaced with the spin cur- 
rent operator, and vice versa. This mapping is useful to 
interpret U < results with the knowledge of the model 
with U >0. 

In a typical TdDMRG implementation the time evo- 
lution operator is represented by the Suzuki- Trotter de- 
composition and the evolution operator at each step is 
efficiently operated to the wave function within an opti- 
mal truncated Hilbert space. Throughout this paper the 
TdDMRG calculations are performed keeping m — 1200 
states and using the 2nd order Suzuki- Trotter decompo- 
sition with time step Ar ~ 0.05h/t. Convergence of the 
results in the limit m —^ oo and At — > is checked for 
U/t ^ 3 and eV/t = 1 (see Fig.l). 

Current between the left and right chains is defined as 
J{t) = e(V'(T)|A^fl|'0(r)) and its time dependence after 
the switching-on of the bias voltage is shown in Fig.l. 
For each parameter set the current starts from 0, ex- 
hibits a certain transient behavior and relaxes to a steady 
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Fig. 1. (Color online) Time dependence of the current after the 
switching-on of the bias voltage for L = 120, Ar = 0.05 and 
m = 1200. Horizontal lines in the figure for U = are the exact 
steady currents for L = oo calculated using Keldysh formalism 
and show that the steady currents can be accurately obtained 
by the TdDMRG calculation. For U/t = 3 and eV/t = 1 we also 
plot the results for (At, m) = (0.025, 1200) and (0.05, 1500) to 
check the two types of errors, the Trotter error and the truncation 
error. The results are almost converged for both Ar and m. 



value accompanied by an oscillatory behavior in some 
cases. After a certain time the steady-like behavior of 
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the current is disturbed. We hereafter call the nonequi- 
librium state in the time interval where the current shows 
a steady-like behavior as a quasi-steady state. 

Termination of the quasi-steady states is caused by the 
reflection of the current at the edge of the system^' ^^' ^^^ 
and thus a finite size effect. Excitations generated by ap- 
plying the bias voltage propagate from left to right with a 
certain velocity and are reflected at the edge, then prop- 
agate back to the left. Eventually those reflected excita- 
tions arrive at the center and disturb the steady flow of 
the current. Because of this effect, for an accurate de- 
termination of the steady current one has to take long 
enough system size to realize complete relaxations before 
the disturbances. However at the same time this effect en- 
ables us to estimate the velocity of the wave front of the 
excitations from the time until the quasi-steady behavior 
ends. 

From Figs.l, we find that V dependence of the ve- 
locity of the wave front is small in the parameter range 
investigated. In the noninteracting case, excitations of 
the system are described by the single particle-hole exci- 
tations in the — 2icosA: band, and the Fermi velocity is 2 
in the units of Fig. 1. Concerning the U > results we do 
not observe significant difference between the velocities 
for eV > Ac and eV < Ac- Thus the excitations which 
carry the initial transient current are almost the same 
for both insulating and quasi-steady states. The veloc- 
ity slightly increases from 2 by increasing U at least for 
U/t < 3. This corresponds to the fact that the velocity of 
the charge excitations of the ground state at half-filling 
is an increasing function of U (see Table I). 

On the other hand for [/ < we see a substantial 
decrease of the velocity. As we stated before, the charge 
excitations of the attractive case correspond to the spin 
excitations of the repulsive case. The observed values in 
Fig.l are 1.5 for U/t = —2 (not shown in the figure) and 
1.3 for U/t = —3. These numbers qualitatively agree with 
the exact values of the spin velocity listed in Table I. 

The oscillation the quasi-steady states exhibit is a fi- 
nite size effect because its amplitude is proportional to 
1/L.^'^°) Note that the effect of the interaction in the 
whole chain strongly influences the amplitudes of the 
oscillations: positive U suppresses the oscillation while 
negative U enhances. This fact suggests that a repulsive 
interaction suppresses coherent transports in the chain 
while an attractive interaction promote them. One in- 
teresting point is that the frequency of the oscillation 
is given by eV for both [/ = and U < 0. For U < 
the initial state has superconducting correlation and thus 
one might expect that the AC Joscphson effect may be 
observed. However the frequency turns out to be always 
eV, not 2eV. The AC Joscphson effect is possible when 
the phase coherence is well developed in each of the two 
subsystems which are weakly linked. In the present sys- 
tem, the phase coherence is difficult to develop because 
the entire system is finite and furthermore the left and 
right subsystems are strongly coupled, t' ~ t. Addition- 
ally for L'^ < and relatively large voltage the current 
shows a damped oscillation. These effects of interaction 
on the oscillatory behaviors provide an interesting sub- 
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Table I. The charge gap Ac, the velocity of the charge excitations 
and the velocity of the spin excitations of the half-filled 1-D 
Hubbard model, calculated from the exact expressions.^^) 
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Fig. 2. (Color online) Current as a function of time for eV/t 
0.2, U/t = 2 and various L, using At = 0.05 and m = 1200. 



ject to be investigated in the future but here we will 
concentrate on the current in the quasi-steady states. 

We show results of J{t) for different system sizes in 
Fig. 2. No essential difference is found in J{t) as long 
as the system retains the quasi-steady behavior. This 
means that the relaxation processes to real nonequilib- 
rium steady states can be well simulated by the finite 
size calculations and that the currents of the steady 
states are obtained from the fiat region. In practice the 
steady currents are estimated from the results of J{t) 
for L ~ 120 inside an interval [30, 50] by the following 
processes: for [/ > taking average, and for t/ < 
fitting data points to the damped oscillation function 
J(r) ~ J{V) + AJe-^/^d'""p sin{eVT + 6). The steady 
currents for [/ > obtained in this way are slightly 
overestimated because the relaxation processes are not 
completely finished in the interval, but its qualitative 
behavior is not influenced by this treatment. 

In order to compare I-V characteristics of the Mott 
insulator with those of band insulator, we also calcu- 
late nonlinear currents in a (noninteracting) band insu- 
lator obtained by replacing the on-site interaction term in 
eq.(2) by an alternating potential (Ab/2) ^ ■ (— l)^njo.. 
This term modifies the dispersion relation from e^ = 
— 2tcosfc to ±-\/e^ + (Ab/2)2 and opens a band gap A;,. 

In Fig. 3 we show the I-V characteristics scaled by the 
gap A, which is Ac for the Mott insulator and A;, for the 
band insulator. By the definition of the energy gap, when 
eV > A the charge excitations are allowed, resulting in 
a finite current. This behavior, dielectric breakdown, is 
beautifully reproduced by our results. Furthermore, in 
each case J{V) for different values of A form a single 
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Fig. 3. (Color online) I-V characteristics of the Mott insulator 
(the 1-D Hubbard model) and the band insulator (see text for 
explanation of the model) scaled by the energy gap A. The re- 
sults of the band insulator are exactly calculated using Keldysh 
formalism and the limit L — > oo is taken. The results of the 
Mott insulator are obtained via the TdDMRG calculation for 
L = 120. Values of the band gap A;, are chosen to be the same 
as the charge gap Ac. (inset) The same I-V characteristics in a 
larger scale. 
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Fig. 4. (Color online) I-V characteristics of the attractive Hub- 
bard model. 



curve when the voltage is smah compared with the band 
width D ~ At. When eV becomes comparable to D the 
energy dependence of the density of states of the band 
is not negligible, and therefore J{V) deviates from the 
scaling curve. For A ^ eV -^ D the scaling curve shows 
almost linear dependence, as shown in the inset of Fig. 3. 
For the band insulator the slope is 2e'^/h, which is the 
value of the perfect conductance. On the contrary for 
the Mott insulator the slope is suppressed by the electron 
correlation and is 1.6e^/h. This surprising scaling behav- 
ior of the 1-D Hubbard model is possible because Ac is 
exponentially small so that all current carrying states are 
under the influence of strong correlation, when eV <^ U. 
We show numerically obtained I-V characteristics for 
U < in Fig. 4. Our results clearly indicate that the lin- 
ear conductance for any t/ < is always the same as the 
perfect conductance 2e^/ft.. The gaplcss charge degrees of 
freedom of the Luther-Emery liquid are described by the 



usual Luttinger liquid theory and the conductance of the 
Luttinger liquid is rcnormahzcd by the correlation expo- 
nent Kp as G = {2e'^/h)Kp.^^^ When we use the Shiba 
transformation eq.(3), Kp in the expression is actually 
given by the correlation exponent for the spin channel 
K^ of the repulsive model. Then if we take the linear 
limit V ^ the model restore the spin SU(2) symmetry 
and this ensures K^r = 1. Accordingly, the linear conduc- 
tance should be equal to 2e^/ft. and this agrees precisely 
with our results. Away from the linear response regime 
J{V) deviates from 2e'^V/h by increasing the voltage and 
the deviation becomes larger with increasing \U\. 

In summary, we have studied nonequilibrium transport 
phenomena in the 1-D Hubbard model at half-filling with 
a finite bias voltage using the TdDMRG technique. We 
have determined the I-V characteristics and found that 
the current for U > in the region eV <^ min(Z?, U) 
shows a universal behavior while the linear conductance 
for [/ < is the perfect conductance 2e^/h. These truly 
nonequilibrium properties were numerically addressed 
with sufficient accuracy for the first time. We believe 
that these reliable data about the nonequilibrium trans- 
port of the 1-D Hubbard model provide basis for future 
studies and stimulate, in particular, analytic approaches. 
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